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Abstract 

We present results from a detailed simulation of a quasi-2D dissipative granular gas, kept in a 
non-condensed steady state via vertical shaking over a rough substrate. This gas shows a weak 
power-law decay in the tails of its Pair Distribution Functions (PDF's), indicating fractality and 
therefore a tendency to form clusters over several size scales. This clustering depends monotonically 
on the dissipation coefficient, and disappears when the sphere-sphere collisions are conservative. 
Clustering is also sensitive to the packing fraction. This gas also displays the standard non- 
equilibrium characteristics of similar systems, including non-Maxwellian velocity distributions. The 
diffusion coefficients are calculated over all the conditions of the simulations, and it is found that 
diluted gases are more diffusive for smaller restitution coefficients. 
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I. INTRODUCTION 



Two-dimensional (2D) granular gases have been widely studied as examples of out-of- 
equilibrium systems that are both simple enough to analyze and easy to construct experi- 
mentally. In the pioneering work of Olafsen and Urbach l|, |2| it was found that these gases, 
kept in a steady state via vertical shaking, could condensate into an hexagonal solid phase 
(for monodisperse sample) as the amplitude or frequency of shaking is reduced. This gives 
a fascinating example of fluid-solid transition not driven by molecular attraction or entropy 
maximization. In the gaseous phase, and confining the vertical expansion of the granular 
layer, the 2D Velocity Distribution Functions (VDFs) displayed by this system are non- 
Maxwellian, and tend to fall instead into an stretched exponential form. The deviation from 
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6[4l2[. and in experiments 



Dehavior of the VDFs in granular gases has been studied theoretically js-S], in 



171], and have been characterized using Sonine 



In its simpler form, experiments in granular gases are carried on using a horizontal cell 
with a perfectly flat bottom; due to dissipation, however, horizontal components of momen- 
tum tend to decrease an thus the gas finally undergoes an inelastic collapse into a static 
(with respect to the plane) condensed phase. In shaken experiments, for large vertical 
acceleration, and in the absence of vertical confinement, there is some dispersion in the 
instantaneous vertical position of the spheres, and the resulting off-plane collisions serve as 
a way of converting vertical momentum into horizontal momentum, keeping in this way the 
gas from condensing. This effect is large enough that it allows for the gas to keep some 
horizontal motion even in cases where a fraction of the beads never loose touch with the 
bottom In recent years some attempts have been made to provide mechanisms by which 
this transfer from vertical to horizontal momentum can be effected without depending on 
the fluctuations in height of the beads. For instance many theoretical and numerical models 
lave been proposed based on an homogeneous granular gas randomly driven by a white noise 
sl, [l8|. Experimentally, a complete layer of spheres glued to the bottom plate (the "floor"), 
or a flrst layer of heavy spheres on which a layer of light spheres is placed, have been used 
, and in other cases an artiflcially roughened bottom plate was employed [l7|, |20|. It 
was found in these experiments that the behavior of the gas was somewhat closer to that of 
an ideal hard disk gas, including velocity distributions closer to the Maxwellian. The use of 
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a roughened substrate makes the movement of the particles to resemble Brownian motion, 
because of the frequent scatterings with the uneven floor. 

Motivated by the experimental results reported in 20|, in this paper we are describ- 
ing a fully three-dimensional simulation of a confined vertically-shaken granular gas, where 
structure has been imposed in the bottom plate using fixed and non-overlapping small hemi- 
spheres. We find that this mechanism is quite efficient in avoiding the collapse of the gas into 
the quiescent state (that is, a state with no horizontal motion). We record the Mean Square 
Displacement (MSD) of the spheres in the gas, and study the dependence of the diffusion 
coefficient on packing fraction, adimensional acceleration and on restitution coefficient. We 
also report the obtained horizontal velocity distributions, and show the strong effects that 
friction with the upper confining plate (the "ceiling") have in these functions. Finally we 
characterize the instantaneous state of the gas using its Pair Distribution Functions (PDF's) 
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23| . which give us information about the underlying effective interparticle interac- 



tions. We observe evidence of week clustering in these PDFs, and relate them to the changes 
in diffusion in the system. 



II. SIMULATION 



We perform numerical simulations of a monodisperse granular gas of spheres in a quasi- 
bidimensional space, square with side L and with periodic boundary conditions in the hori- 
zontal. This gas is confined between two horizontal planes which oscillate in the vertical in 
a sinusoidal way z = A sin(a;t), with amplitude A and frequency uj large enough so that the 
spheres will in general touch both confining plates in every oscillation of the system. In the 
bottom plate there are fixed and non-overlapping hemispheres whose diameter cx/js is almost 
half of the diameter of the free spheres as- These hemispheres move synchronously with the 
plane and are randomly distributed with a sufficiently high two-dimensional packing fraction 
(phs, so that no wide fiat patches can be formed (see Fig. ([1])). Here (phs = ^hs7rcr^g/(4L^), 
where N^g is the number of hemispheres in the plate of side L. The separation between 
planes is h = l.Gag + (Jhs/4: ~ I-TcTs, and the minimum distance between the centers of two 
hemispheres is dm- 

In order to solve the Newton's equations of motion we use standard time-driven molecular 
dynamics, using a velocity- Verlet integration with a predictor step for velocities. This is nec- 
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essary in order to incorporate the dissipative part of the force, which is velocity-dependent. 
The force between particles in contact is described by 

F,, = I + ^ °' (1) 

I otherwise, 

where F^- is a normal force which causes changes of the translational motion of the particles 
and Fjj is a tangential force, originating in friction, which causes changes in both the trans- 
lational and the rotational motion. The quantity is the mutual overlap (compression) of 
particles i and j and is defined by 

^ij = max(0, as - \ri - rj\) (2) 

for interparticle collisions, by 

^ij = max(0, {as + ahs)/2 - |ri - Rj|) (3) 

for particle-scatterer collisions, and by 

^ip = max(0, as/2- dip) (4) 

for sphere-plate collisions. Here denotes the center of the i-th grain, Rj the center of the 
i-th scatterer, and dip is the normal distance from the center of the i-th grain to the surface 
of the plate p. 

Modeling a force that leads to inelastic collisions requires at least two terms: repulsion 
and some type of dissipation. The existiiig models and their characteristics, as well as 



comparisons among them, can be found in 2J-[29|. Here we have used for the normal force 



the linear spring- dashpot model, in which the contact interaction is modeled by the damped 
harmonic oscillator force 

F„ = min jo, - 7„^| , (5) 

where 7„ is a damping constant and k„ is related to the stiffness of a spring whose elongations 
is ^, the overlap between two grains (or a grain and a boundary plate). The advantage of 
this model lays in its analytic solution, where the collision time is given by 

] 

2mefr, 




tn = -{i^-{T^] 1 , (6) 



4 



and the restitution coefficient is 

Cn = exp ( - tn] . (7) 

Here is the effective mass for the colhding pair. For the tangential force we consider 
the Coulomb dynamical friction law F-^- = — /iF^"v*^-, where vjj is the unitary vector in the 
direction of the relative tangential velocity between spheres i and j (or, by extension, between 
a sphere and a scatterer or a sphere and a plate). The effect of gravity is incorporated into 
the vertical acceleration. 

For the actual simulations we have taken most parameters similar to those of the pre- 



viously mentioned experiments 20|], as performed with steel spheres. Therefore, we have 
fixed the frequency of oscillation of the cell to / = 60 Hz, and have considered two ampli- 
tudes of oscillation: Ai = 0.024 cm and A2 = 0.05 cm, for every packing fraction of free 
spheres, getting in this way for the adimensional acceleration F = (2'KfYA/g the values 
Fi = 3.5 and F2 = 7.2. Here is defined as before by = A^crf 7r/(4L^). As parameters 
for the normal force, Eq. ([5]), we use those extracted from Eqs. fl6ll71) . and fix the collision 
time to tn = 7.1x10"^ s, and use three values of the restitution coefficients e„ = 0.36,0.66 
and 1 for the sphere-sphere and sphere-hemisphere interactions (at each value of F and 
for every packing fraction considered). The coefficient of restitution for sphere-substrate 
interaction was fixed at 0.878, corresponding to an experimentally measured steel-acrylic 



restitution coefficient 20], and the considered coefficient of dynamical friction fi was 0.25. 
The time-step of the integration is always fixed to t„/100, and the simulation time is 40 s, 
after a transient of 10 s. The diameter of the free spheres ag and of the hemispheres ahs have 
been fixed to 0.44 cm and 0.2 cm, respectively, in all simulations. The minimal distance 
between hemispheres is = 0.2 cm. Finally, the boundary plates have side L = 16 cm, and 
the packing fractions used have been = 0.15, 0.20, 0.25, 0.30 and 0.35 for both F values 
used. The resulting numbers of grains are Ng = 252, 336, 420, 504 and 588 respectively, and 
the number of scatterers is N^s = 4100. All the quantities reported here are the average 
over 10 simulations, each with a different configuration of hemispheres and different initial 
conditions. 
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III. NUMERICAL RESULT: DIFFUSION 



We will start by presenting results for the diffusion in the quasi-2D granular gas. As 
mentioned before, the interaction with the hemispheres provides horizontal momentum to 
the particles, but also scatters their motion, making them to move on a pseudorandom way, 
resembling in certain form a Brownian motion. It is clear that the interparticle collisions also 
contribute to the scattering. In Fig. (|2]) we can see two trajectories for a single particle on 
a complete simulation of 40 s, one for the case with e„ = 1 and another for e„ = 0.36. This 
was done with a packing fraction cj) = 0.15, and the trajectory includes the interactions with 
the other particles and with the substrate. One must remember here that the restitution 
coefficient for particle-plate collisions is always fixed, while that for particle-hemisphere 
collisions is the same as the one for particle-particle collisions, and can be varied. 

The MSDs measured fit very well the Einstein form 

((r(t + to)-r(to))')=4Dt; (8) 

this is not as trivially expected as it may seem, since some recent results point to a breakdown 
of Einstein's law for 2D granular gases [s^. This is related to a well-known anomaly for 
self-diffusion in 2D, where the presence of hydrodynamical backfiows gives origin to long 
tails in the velocity autocorrelation function. These tails behave as t~^, giving in this way a 
ogarit 
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imic divergence to the diffusion coefficient, accordingly to the Green-Kubo formula 
In our simulations we have not seen any evidence of non-Einstenian behavior; in 
particular we have performed sets of longer runs (80 s) for three different sets of parameters, 
and one set of runs for a bigger system (L = 16a/2 = 22.63 cm), and in all cases the same 
diffusive behavior was found, with the diffusion coefficient independent of the length of the 



run or the system size. Looking at the difference between these results and those of [30|, it 
is clear that the presence of fixed (with respect to the plane) scatterers is the reason why 
regular diffusion is restored. 

In Fig. ([3]) we show the diffusion coefficient D versus the packing fraction for both values 
of the adimensional acceleration. Here we can observe the following points: first, D increases 
for low values of (t>, and decrease for high values of 0. These results are in agreement with 



those obtained by [33|, 13^ , and are to be expected, since larger densities represent in general 
more obstacles to the motion. Second, for small values of the diffusion coefficient increases 
when we lower e„, and this behavior is stronger for low F. As the concentration and F grow, 



this behavior reverses, that is, for large there is a crossover were the diffusion grows 
with e„. Theoretical evidence for an inverse relationship between diffusion and restitution 
coefficients has recently been found in [33], as shown in their Fig. (5). This is also noticeable 
in our Fig. ([2]), where one can observe how the particle excursions are larger for lower e„. 

The increase in diffusion for smaller e„ is related to the appearance of density fluctuations 
that become stronger as the dissipation increases. This weak clusterization becomes apparent 
as a fall of the tails of the PDFs for large distances, as we will see later. Clusterization 
ends up freeing some space in the system, allowing in this way faster diffusion. The same 
phenomenon is well known in colloid-polymer mixtures, where the depletion forces induced 
by the (small) polymers create strong density fluctuations in the colloidal phase, and increase 



its diffusion 
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IV. GRANULAR TEMPERATURE AND VELOCITY DISTRIBUTION FUNC- 
TIONS 

In this section we perform a detailed analysis of the granular temperature Tg and of the 
Velocity Distribution Functions (VDFs) P{v) obtained in the simulations. In Fig. (jl]) we 
can observe the behavior of the in-plane Tg, defined as Tg = T^ + Ty = {vl) + {Vy), as a 
function of the packing fraction 0, for all values of F and e„ considered. The following points 
deserve to be mentioned: first, for any given values of and^, Tg increases together with 



F, as is intuitively expected, and experimentally observed [17|. Second, one can observe 
that for given values of and F, Tg decreases together with e„, also as expected. Third, 
and more interesting, we find that for^ = 1 the Tg grows slowly with 0, in agreement with 



the experimental results obtained in 17|]. In that experiment Tg grows monotonously until 
reaching a maximum for some 0max ~ 0.5, and then decays as keeps growing. We do not 
find this later decay of Tg, probably because we are using values of below 0max- 

The coincidence found for this particular value of e„ deserves some discussion: although 

n fl n 

the authors of [17] do not report their working value of e„, other references |37[ give 
Cn ~ 0.9 for the stainless steel spheres used there. As relates to the simulation, we remind 
the reader that even when e„ = 1 there is some dissipation in our system, since this value of 
the restitution coefficient is applied only for sphere-sphere and sphere-hemisphere collisions. 
Contacts with the upper plate (and, even if infrequent, with the lower plate), are still 
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inelastic, with = 0.878. Besides, the existence of a nonzero value of /i gives us some 
energy losses via friction, for all collisions. One may therefore assume that there is, in our 
simulation, some effective e„ a bit smaller than 1, and that the slow growth of the Tg for 
diluted gases found in both experiment and simulation occurs for high values of e„, not 
necessarily e„ = 1. This should be stressed since, for values of e„ well bellow 1, Tg actually 
decays very slowly as increases, at least for the range of considered here. We are not 
aware of any experimental study of Tg for highly dissipative materials, and so there is no 
verification of this curious behavior reversal for small restitution coefficients. 

The VDFs for both components of the horizontal velocity have been obtained, for all 
values of F, and e„ considered. In Fig. we show a few of those curves, normalized 
by their characteristic velocities Vc = \/ (vf) = y/Ti. In the same Figure we also show, 
as a reference, a unitary Gaussian. The first remarkable fact is that for all values of 0, 
when e„ = 0.66 and 1.0, with F = 3.5, and also for all values of e„ and F = 7.2, the 
VDFs present a strong peak at the center of the distribution — that is, for low velocities — . 
(Only some examples are shown, to avoid crowding the Figure.) These peaks can be also 
found in Fig. (8) of jl^, for small 0. In that reference, they are attributed to the fact 
that at these packing fractions the effect of interparticle collisions is minimal, and so the 
general distribution is dominated by the VDF corresponding to one isolated particle. The 
simulations show that the main reason for the appearance of this peak is the friction of the 
grains with the upper plate in the system. In particular, as can be easily visualized in a 
ID system, rotational velocities acquired at the collisions with the floor act as a braking 
factor upon contact with the ceiling, increasing in this way the concentration of particles 
with low horizontal velocities. In Fig. ([6]) we show the normalized curves for the VDFs, for 
the case = 0.20, F = 3.5, and e„ = 0.66 and 1.0, but taking the grain-ceiling value of fi 
as 0.0 in one case and /i = 0.25 in the other. It is quite apparent from this Figure how, 
when there is no friction with the ceiling, there are no low-velocity peaks. However, it can 
be seen from the same Figure that the braking effect of the ceiling is not the only factor 
involved in the non-Gaussian behavior of the VDFs; in particular, the line for e„ = 0.66 
and /iceiiing = shows very clear deviations from the Gaussian. Even stronger deviations 
are found for e„ = 0.36 and /iceiUng = (not shown in the Figure.) We would also like to 
mention that these peaks diminish as the packing fraction of the systems grows, as observed 
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To finish this section, let us show some analysis of the tails of the VDFs. In Figs. ([7]) 
and ([8]) we plot — log(— log(P(f /fc)/P(0))) vs. log(f/fc) for F = 3.5 and e„ = 0.36, and 
for F = 7.2 and e„ = 1.0, as a function of 0. For a perfect Gaussian centered at the origin, 
this plot should give a straight line with slope m = —2. It can be observed from both 
Figures how the behavior of the tails of the VDFs change slowly with 0, but quite strongly 
as IS varied. We include in these graphs linear fits for the intermediate and final parts 
of the distribution. Here it is important to notice that for F = 3.5 y e„ = 0.36, both the 
intermediate and final parts of the VDFs are close to a Gaussian, as can be seen also, for 
(p = 0.15, in Fig. However, when e„ = 1.0, we find the intermediate and final parts of 
the distribution to be clearly different. While the final part is reasonably close to a Gaussian 
[nit ~ 2) the intermediate part has an exponent larger than —1. This is due to the effect of 
the ceiling, since the deformation of the VDFs — that is, their deviation from a Gaussian — 
in this region is related to the increase in the number of slow particles, which is in turn due 
to the friction with the upper plate. This general behavior can be corroborated in Fig. ([8]) 
(F = 7.2, Cn = 1.0). Here we can see how the final part of the distribution remains lose to 
Gaussian, while one gets an exponent rrit ~ —0.5, much larger than —2, for the intermediate 
part of the distribution. 

One then finds that the effect of the rough substrate is to return the form of the VDFs to 
the Gaussian form, since the presence of the scatterers breaks translational symmetry and 
also forces a better horizontal energy interchange between the particles. On the contrary, 
the effect of the friction with the ceiling creates an increase in the distribution for slow 
particles, with a corresponding drop for the distribution at intermediate velocities. 



V. PAIR DISTRIBUTION FUNCTIONS 



The Pair Distribution Functions (PDFs) for a packing fraction = 0.35 are shown 
in Fig. (Q, for some combinations of F and e^. It can be observed how when e„ = 1, 
the PDFs do not depend on F, and the structure observed is very similar to that of an 
equilibrium elastic hard sphere gas j2|, indicating that the correlations that exist are mostly 
due to excluded volume effects. On the other hand, when the dissipation between particles 
increases, the correlations between particles begin to grow, an effect that is more notorious 
for low values of F. 
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The lines shown in the graph have two intriguing characteristics: first, we find a complete 
absence of any structure at r/cr = ^/S = 1.73, where a peak appears in the close-packed limit. 
Neither there is a peak around r/a = \/2 = 1.41, which would signal some type of square 
clustering. It is therefore clear that the rough floor eliminates completely the possibility of 
formation of large crystalline clusters. Second, a very peculiar secondary peak develops, for 
small driving and high dissipation (T = 3.5 and e„ = 0.36), at r/a ~ 2. The possible reason 
we have identifled for the origin of this peak is the formation, at this low values of driving 
and restitution coefficient, of short-lived linear chains of grains (see Fig. ([1]) for an example). 
These linear chains may be induced by some residual order in the substrate, an order that 
in turn is due to the fact that it has to satisfy simultaneously a large density of scatterers 
and the no-overlap condition. It will be therefore interesting to see the evolution of q(r) 

nn 

in those experiments where an ordered substrate has been used [15|, 119|; although it is also 
clear from the graph that only for the smallest restitution coefficient used we get this peak. 

Fore a more comprehensive view of the evolution of g{r), in Fig. ( ITOj) we show the PDFs 
for r = 3.5, as a function of e„ and 0. The very strong effect that e„ has in these curves is 
the first thing to notice; for high dissipation the structure of the main peak in g{r) is almost 
independent of (p, and little variation can be seen in the rest of the curve. The large value of 
g{r) at the main peak signals a tendency to clustering. At the other extreme, the gas with 
e„ = 1.0 shows much more sensitivity to (p. An unusual behavior can be seen for e„ = 1 
and small 0: the largest value of g{r) no longer falls at contact, and instead, a soft "hill" 
develops around r/a ^ 1.2-1.3. This is again an effect of the scatterers, that a this high 
value of e„ intrude between the grains, and, for small packing factors, reduce the probability 
of contact. Again, there is no signal of any peak at r/a = v^, but we find instead that for 
the most dissipative system, some peak at r/a = 2 begins to form. Visual inspection of the 
dynamics in this parameter sector shows in effect the presence of short-lived linear chains. 
For r = 7.2 (not shown) the results found are analogous: there is an increase in the height 
of the main peak of g{r) clS IS lowered, although not as strong as for F = 3.5; and there 
is absolutely no peak at r/a = As for the behavior around r/a = 2, now only some 
weak increase in g{r) can be detected. A remarkable result is the fact that for all (p used 
the different g{r) obtained for e„ = 1 agree completely for both values of F. 

Finally, one of the most interesting points in the phenomenology of this granular gas is 
shown in Fig. f llip . which displays a close up to the PDFs for F = 3.5. It can be observed 
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here that for each filhng fraction considered, g{r) has its long-r tail falhng shghtly below 
1.0 when we lower the restitution coefficient, something that does not happen in the curves 
with e„ = 1. We would like to stress that the same normalization protocols were used for all 
lines. This behavior implies a weak fractal behavior related to the non-unitary restitution 
coefficients in the system. More in detail, one can observe that the tails in the PDFs show 
the largest drop below 1.0 for = 0.15 and e„ = 0.36; in general, the PDFs approach 1.0 
from below as (p increases or increase. As an effect of the clusterization signaled by the 
falling of the PDFs bellow 1, there is an increase in the mean free path of those particles 
outside of the clusters — that is, in the more diluted part of the gas — . This gives rise to 
an increase of the average diffusion coefficient of the system, as was shown in Fig. For 
r = 7.2 one gets some appreciable decay in the tails of g{r) only for small (p and small e„, 
and in general this effect is not so strong. This behavior again coincides with the increase 
in diffusion found for this value of F, as can be seen in Fig. (|3]). 

It is important to notice that the decay observed in the tails of the PFDs is due only to 
interparticle collisions, and independent of friction, either with the floor or the ceiling. This 
can be observed in Fig. ( ITTj) . where the fastest decay of the tails of the g{r) appears for 
r = 3.5 y e„ = 0.36; for these values of the parameters the VDFs are practically Gaussian, 
as can be seen in Fig. (|5]), implying very little interaction with the ceiling. 

So, it is clear, looking at the results for diffusion and for the PDFs, that clustering 
increases the average diffusion, even if these clusters are very unstable and their lifetimes 
are very short. The weakly fractal behavior of the gas disappears when we increase any 
of the three controlling parameters in the system: the packing fraction, the amplitude of 
vibration or the restitution coefficient. 

VI. CONCLUSIONS AND FINAL COMMENTS 

We have performed numerical simulations of a vibrating quasi 2D granular gas, using 
a fully 3D algorithm with Coulomb friction proportional to normal forces, for a horizontal 
granular gas confined by a rough floor and a flat ceiling. The system shows some of the 
characteristics already explored in experimental and other simulational work, as for instance 
the recovery of Gaussian VDFs at low F over a rough substrate. The principal result of the 
simulation is how, even with a floor full of scatterers, there is some level of clustering for 
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the low range of F, and values. This clustering shows in the strong growth of the main 
peak of g{r), in the increase in diffusivity lowered, and very clearly in the drop 

of g{r) below 1 for for r ^ a. All of the characteristics that indicate clustering coexist, 
however, with an almost complete recovery of Gaussian behavior in the VDFs for this part 
of the parameter space. For the other corner of this space — that is, high F and and 
e„ = 1 — all signals of clustering disappear, making it clear that the rough floor is only of 
minor importance with respect to this behavior. 

There are also some very punctual phenomena associated to this simulation: first, it 
shows quite dramatically the distorting effect that the ceiling may have on the VDFs, by 
introducing a braking affect that results in a large increase in the central part of the dis- 
tributions — the very low velocity regime — . Experimental results reported up to now show 
some hints of this problem, but there is not, to our knowledge, any systematic experimental 
exploration of this effect. Also, the presence of a peak in g{r) around r /a — 2 is pointing 
towards some complex effects of the substrate on the dynamics. Even if the configurations 
of scatterers used in the simulation have a random origin, the simultaneous demands large 
density and no overlap may have introduced some level of ordering, later reflected is a 
larger-than-expcctcd population of linear and almost-linear chains in the granular gas. The 
appearance of undesired order in the substrate is not surprising, since it is actually quite 
difficult to pack monodisperse disks densely without generating a crystal. 

Acknowledgments 

This work was supported by CONACyT (Mexico) through Grant No. 82975. J. A. P.-B. 
receives support from a GONAGyT Doctorate Fellowship. 



[1] J. S. Olafsen and J. S. Urbach, Phys. Rev. Lett. 81, 4369 (1998). 
[2] J. S. Olafsen and J. S. Urbach, Phys. Rev. E 60, R2468 (1999). 

[3] A. Puglisi, V. Loreto, U. Marini Bettolo Marconi and A. Vulpiani, Phys. Rev. E 50, 5582 
(1999). 

[4] N. V. Brilliantov and T. Poschel, Phys. Rev. E 61, 5573 (2000). 
[5] T. van Noije and M. Ernst, Granular Matter 1, 57 (1998). 

12 



[6] Y. H. Taguchi and H. Takayasu, Europhys. Lett. 30, 499 (1995). 

[7] X. Nie, E. Ben-Nairn and S. Y. Chen, Europhys. Lett. 51, 679 (2000). 

[8] A. Barrat and E. Trizac, Phys. Rev. E 66, 051303 (2002). 

[9] J. S. van Zon and F. C. MacKintosh, Phys. Rev. Lett. 93, 038001 (2004). 
[10] J. S. van Zon and F. C. MacKintosh, Phys. Rev. E 72, 051301 (2005). 
[11] D. J. Bray, M. R. Swift and P. J. King, Phys. Rev. E 75, 062301 (2007). 
[12] A. Burdeau and P. Viot, Phys. Rev. E 79, 061306 (2009). 

[13] W. Losert, D. G. W. Cooper, J. Delour, A. KudroUi and J. P. Gollub, CHAOS 9, 682 (1999). 
[14] A. KudrolU and J. Henri, Phys. Rev. E 62, 1489 (2000). 

[15] A. Prevost, D. A. Egolf and J. S. Urbach, Phys. Rev. Lett. 89, 084301 (2002). 

[16] K. Kohlsted, A. Snezhko, M. V. Sapozhnikov, I. S. Aranson, J. S. Olafsen and E. Ben-Nairn, 

Phys. Rev. Lett. 95, 068001 (2005). 
[17] P. M. Reis, R. A. Ingale and M. D. Shattuck, Phys. Rev. E 75, 051311 (2007). 
[18] D. R. M. Wilhams and F. C. MacKintosh, Phys. Rev. E 54, R9 (1996). 
[19] G. W. Baxter and J. S. Olafsen, Nature 425, 680 (2003). 

[20] R. A. Bordallo-Favela, A. Rami'rez-Safto, C. A. Pacheco-Mohna, J. A. Perera-Burgos, Y. 
Nahmad-Mohnari and G. Perez- Angel, Eur. Phys. J. E 28, 395 (2009). 

[21] B. V. R. Tata, P. V. Rajamani, J. Chakrabarti, A. Nikolov and D. T. Wasan, Phys. Rev. Lett. 
84, 3626 (2000). 

[22] P. Eshuis, K. van der Weele, D. van der Meer and D. Lohse, Phys. Rev. Lett. 95, 258001 
(2005). 

[23] F. Pacheco- Vazquez, G. A. Caballero-Robledo and J. C. Ruiz-Siiarez, Phys. Rev. Lett. 102, 

170601 (2009). 

[24] J. Schafer, S. Dippel and D. E. Wolf, J. Phys. L France 6, 5 (1996). 

[25] N. V. Brilliantov, F. Spahn, J.-M. Hertzsch and T. Poschel, Phys. Rev. E 53, 5382 (1996). 
[26] A. Di Renzo and F. P. Di Maio, Chemical Engineering Science 59, 525 (2004). 
[27] A. B. Stevens and C. M. Hrenya, Powder Technology 154, 99 (2005). 

[28] H. Kruggel-Emden, E. Simsek, S. Rickelt, S. Wirtz and V. Scherer, Powder Technology 171, 
157 (2007). 

[29] T. Poschel and T. Schwager, Computational Granular Dynamics: Models and Algorithms 
(Springer, 2004). 

13 



[30] C. Henrique, G. Batrouni and D. Bideau, in Granular Gases, edited by T. Poschel and S. 

Lading, Lecture Notes in Physics vol. 564, p. 140. 

[31] B. J. Alder and T. E. Wainwright, Phys. Rev. A 1 18 (1970). 

[32] P. J. Camp, Phys. Rev. E 71 031507 (2005). 

[33] W. T. Kranz, M. Sperl and A. Zippelius. [a"rXiv:1002.2150V l. 

[34] A. R. Abate and D. J. Durian, Phys. Rev. E 74 031308 (2006). 

[35] K. N. Pham et al, Science 296 104 (2002). 

[36] R. Juarez-Maldonado and M. Medina-Noyola, Phys. Rev. E 77 051503 (2008). 

[37] A. Kudrolh, M. Wolpert and J. P. Gollub Phys. Rev. lett. 78 1383 (1997). 



14 




FIG. 1: (a) Fraction of a snapshot (taken from above) of the granular gas, covering around 1/13 of 
the actually simulated area. The hemispheres (scatterers) are shown as white circles, the moving 
grains as yellow circles. Notice the presence of some short linear chains, (b) Schematic lateral view 
of the simulated system, to scale. 
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FIG. 2: Trajectory of a single particle for two different restitution coefficients. The more compact 
trace (blue) is for e„ = 1.0, the less compact one (green) is for e„ = 0.36. The packing fraction is 
0.15 and we used T = 3.5. Notice how the tracer describes a much larger excursion for the smaller 
restitution coefficient. 



16 




FIG. 3: Diffusion Coefficients vs Filling Fraction for: (a) F = 3.5, and (b) F = 7.2. The values of 
e„ are, at the left and starting from top: e„ = 0.36 (red), = 0.66 (green) and = 1 (blue). For 
the range of packing fractions covered, D always decreases with (/). The behavior with respect to 
e„ is more complex: for small (p, D decays as e„ grows. For large (f) there seems to be a crossover 
where D starts growing with. Cn. 
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FIG. 4: Tg vs (j) for both values of (a) T = 3.5 and (b) T = 7.2. The values of e„ are, for both 
graphs and starting from top: e„ = 1 (blue), e„ = 0.66 (green), e„ = 0.33 (red). The dependence 
of Tg on (p is much weaker than its dependence on Cji. For both values of T we find Tg decaying 
slowly with cf) at low values of e„, but increasing, agin slowly, with </> for e„ = 1. 
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FIG. 5: Velocity Distribution Functions (VDFs) for some values of the parameters. Here the VDFs 
have been scaled by the RMS velocity. Looking at the center, and starting from the top, the four 
thick lines correspond to: F = 7.2, cj) = 0.15 and e„ = 0.36 (magenta); F = 3.5, = 0.15 and 
e„ = 1.0 (green); F = 3.5, = 0.35 and e„ = 1.0 (blue); and F = 3.5, 4> = 0-15 and e„ = 0.36 
(red). The thin black line corresponds to a Gaussian distribution. For the last set of parameters 
the combination of high dissipation and small F results in trajectories where the grains almost 
never touch the ceiling, giving therefore a VDF very close to the Gaussian. 
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FIG. 6: Some VDFs with and without friction with the ceihng. For all lines we have used here 
r = 3.5 and (p = 0.2. Looking at the center, and starting from the top, the four thick lines 
correspond to (a) = 1.0 and n = 0.25; (a) e„ = 0.66 and = 0.25; (a) e„ = 0.66 and /i = 0; and 
(a) = 1.0 and /i = 0. The thin black line is again a Gaussian distribution. 
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FIG. 7: Intermediate and outer parts of the VDFs for T = 3.5 and e = 0.36, for different values of 
(j). Starting from top, the values of (j) are 0.15, 0.20, 0.25, 0.30, 0.35. For clarity, lines have been 
shifted apart from each other. Again, low F and Cn reduce contacts with the ceiling, therefore 
generating VDFs very close to Gaussians. The slope m in the intermediate region is m = —1.85 
(thin blue segment), and for the outer region (the tails) we get m = —1.97 (thin black segment). 
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FIG. 8: Intermediate and outer parts of the VDFs for F = 7.2 and e = 1.00, for different values 
of (j). The values of (j) are as in the previous Figure, and again the lines have been shifted apart 
from each other. Here the high F and tn imply strong bounces in the ceiling, and therefore a 
very noticeable shift of the VDF towards low velocities. These VDFs are quite distorted, and 
the double-log vs log graph cannot be fitted well with straight lines. The best linear fits for the 
intermediate and outer regions give m = —0.66 (thin blue segment) and m = —1.94 (thin black 
segment), respectively. 
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FIG. 9: Pair distribution functions for = 0.35 for some values of F and e„. At the main peak, 
and starting from top, the changing parameters are: (a) F = 3.5, e„ = 0.36; (b) F = 7.2, e„ = 0.36; 
(c) F = 3.5, en = 0.66; (d) F = 7.2, e„ = 0.66; (e) F = 3.5, e„ = 1.0 and (f) F = 7.2, en = 1.0. 
The last two lines fall on top of each other, showing how the PDFs for non-dissipative grains are 
independent of the acceleration. In general, smaller restitution coefficients give rise to larger values 
of g{r) at contact. 



23 



1 r 

a) 




1 1 1 

1 \ 

b) 


1 


_ 


_ 






- 


- 


1 ^^^^^^ 

j 

J\ 1 1 


1 




1 1.5 2 2.5 3 

r/a 

FIG. 10: Pair distribution functions for Fi = 3.5, for: (a) e„ = 1.0, (b) e„ = 0.66, and (c) 
Cn = 0.36. In all cases, the g{r) with the largest variation corresponds to <j) = 0.35, and g{r) 
becomes progressively flatter as (j) decreases. Notice that in (a), for (p < 0.30, the largest peak is 
not located at contact, indicating that the fixed scatterers have a weak separating effect over the 
grains. 
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FIG. 11: Zoom of the g{r) for Fi = 3.5. The decreasing of the g{r) tails can be observed for 
different filhng fractions with respect to the curve with e = 1. The parameters are, starting from 
below: (a) = 0.36, (/> = 0.15 (red), (b) e„ = 0.66, (/> = 0.15 (blue), (c) e„ = 0.36, </> = 0.25 
(green), and (d) e„ = 0.66, (p = 0.25 (magenta). The two sets of points running at g(r) = 1 
correspond to e„ = 1, for (p = 0.15 and 0.25. The lines are fits of the form y{x) = A + B/{x — xq), 
and are provided only as guides to the eye; at this reduced vertical scale, the large noise levels in 
the data make a detailed evaluation of the parameters of the fit too uncertain as to be useful. 




25 



